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We investigate the nonequilibrium phenomena through the quantum dot coupled to the normal 
and superconducting leads using a weak-coupling continuous-time Monte Carlo method. Calculating 
the time evolution of particle number, double occupancy, and pairing correlation at the quantum 
dot, we discuss how the system approaches the steady state. We also deduce the steady current 
through the quantum dot beyond the linear response region. It is clarified that the interaction 
decreases the current in the Kondo-singlet dominant region. On the other hand, when the quantum 
dot is tightly coupled to the superconducting lead, the current is increased by the introduction of 
the Coulomb interaction, which originates from the competition between the Kondo and proximity 
effects. Transient currents induced by the interaction quench are also addressed. 

PACS numbers: Valid PACS appear here 



I. INTRODUCTION 

Recently electron transport through nanofabrications 
has attracted much interest. One of the simplest systems 
is a quantum dot with discrete energy levels^ which gives 
us a stage to study fundamental quantum physics. When 
the quantum dot is contacted with the normal leads, 
electron correlations play a crucial role in understand- 
ing its transport at low temperatures where the Coulomb 
blockade and Kondo effect appear. On the other hand, 
when superconducting leads are connected to the quan- 
tum dot^ the proximity-induced on-dot pairing becomes 
important, in addition to electron correlations. However, 
multiple Andreev reflections should dominate the system 
and it is difficult to study the interplay between the su- 
perconducting and electron correlations systematically. 

The quantum dot system coupled to the normal and 
superconducting leads is one of the appropriate systems 
to study how the transport properties are affected by the 
competition between the Kondo and proximity-induced 
on-dot pairing effects. In fact, the system has exper- 
imentally been examined^— and the Kondo-enhanced 
Andreev transport has recently been observed in the 
InAs quantum dot.— i Theoretical study has been done by 
many groups,— ~— and some interesting transport proper- 
ties have successfully been explained. However, it is non- 
trivial how the techniques are applicable in the strong 
coupling and high voltage region. This may be impor- 
tant to understand the experimental results quantita- 
tively since the linear response region is narrow in the 
quantum dot system. Therefore, the unbiased and robust 
method for the nonequilibrium phenomena is desired. 

To this end, we make use of the continuous-time quan- 
tum Monte Carlo (CTQMC) method^ based on the 
Keldysh formalism . 18 ' 19 Here we extend the CTQMC 
method in the continuous-time auxiliary field (CTAUX) 
formulation^ to treat the superconducting state in the 
Nambu formalism. Calculating the particle number, dou- 
ble occupancy, pairing correlations and current through 
the quantum dot, we study the nonequilibrium phenom- 



ena. We also discuss the competition between the Kondo 
and proximity effects on the steady-state in the quantum 
dot coupled to the normal and superconducting leads. 

The paper is organized as follows. In Sec. |Hl we intro- 
duce the model Hamiltonian for the quantum dot coupled 
to the normal and superconducting leads. The CTQMC 
algorithm in the Nambu formalism is explained in Sec. 
TTT1 In Sec. IIV| we discuss the nonequilibrium phenomena 
in the quantum dot system. A brief summary is given in 
Sec. H 



II. MODEL 

We consider the electron transport through the quan- 
tum dot coupled to the normal and superconducting 
leads, which are labeled by a — N and S. For simplic- 
ity, we use a single level quantum dot with the Coulomb 
interaction. The system should be described by the fol- 
lowing Anderson impurity Hamiltonian as 
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where Cfc QCT (cj. Q(T ) is the annihilation (creation) operator 
of an electron with wave vector k and spin cr(=t, 4-) in 
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the ath lead, d IJ {d) IJ ) is the annihilation (creation) oper- 
ator of an electron at the quantum dot and n a = S a d a . 
ek a is the dispersion relation of the ath lead and Vk a is 
the hybridization between the ath lead and the quantum 
dot. ed and U is the energy level and the Coulomb in- 
teraction at the quantum dot. To discuss the nonequilib- 
rium state in the system, we set the chemical potential 
in each lead as //jv = V and us = 0, where V is the 
bias voltage. In our paper, we focus on the particle-hole 
symmetric system with + U/2 = and the supercon- 
ducting lead is assumed to be described by the BCS the- 
ory with an isotropic gap As = A(> 0). We consider a 
sufficiently wide bandwidth in both leads, where the cou- 



pling strength T a (u) = i^Y^k 
constant. 
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When no bias voltage is applied to the system, ground 
state properties depend on the ratio Ts/^n- In the case 
Ts/rAr <C 1 and U ^ 0, conduction electrons in the nor- 
mal lead screen the local spin at the quantum dot and the 
Kondo-singlet dominant state is realized. On the other 
hand, the singlet Cooper pairs are realized at a quantum 
dot due to the proximity effect when r^/IV 3> 1- It is 
known that when the ratio is changed, the crossover, in 
general, occurs between these two singlet states and the 
first-order transition occurs in the limit IV = 0. 14 ' 15 ' 21 



This crossover affects noncquilibrium properties in 
the quantum dot system. It has been reported how 
the zero bias conduct anc o 1 1 1 1 3 and the current-voltage 
characteristics^ are affected in the vicinity of the 
crossover. Furthermore, a detailed structure in the dif- 
ferential conductance has been discussed in the nonlin- 
ear response region by means of the modified perturba- 
tion theory (MPT)^ which is based on an interpolation 
scheme between the weak-coupling limit (U — > 0) and 
the superconducting atomic limit (Tjv = 0, A — > oo). 
Although the reliable results have been obtained in the 
cases U/Tn <C I and hs/r^y 3> I, it is unclear whether 
the nonlinear response in the strong coupling region is 
quantitatively described or not. 



In this paper, we make use of the CTQMC method. In 
the method, Monte Carlo samplings of collections of dia- 
grams arc performed in continuous time and thereby the 
Trotter error, which originates from the Suzuki- Trotter 
decomposition, is avoided. This method has successfully 
been applied to many equilibrium systems. Recently, the 
CTQMC method based on the Keldysh formalism has 
been formulated, where the nonequilibrium phenomena 
in the quantum dot coupled to normal leads have quan- 
titatively been studiedi 18 ' 19 In the following, we extend 
the CTQMC method to deal with the superconducting 
state in the Nambu formalism. 



III. CONTINUOUS-TIME QUANTUM MONTE 
CARLO SIMULATIONS IN THE NAMBU 
FORMALISM 

In this section, we explain the CTQMC method based 
on the Keldysh formalism) 18 i 19 and extend it to treat 
the superconducting state. Here, we consider a weak- 
coupling version of the CTQMC approach. Since the 
interaction is considered as a perturbation, we can exam- 
ine the time evolution of the system after the interaction 
quench. We start with the following identity 
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where po(= e~P H ° /Tre~^ H °) is the initial density matrix 
for Ho(= H — H') and K is some nonzero constant. It is 
then given by 
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where 0{t) = e ltH °Oe ltH ° is the interaction represen- 
tation of the operator O and T(T) is the time-ordering 
(antitime-ordering) operator. Expanding the exponen- 
tials into a power series, we obtain 
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where we have used the following equation as 
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with 7 = cosh _1 (f + tU/2K). The introduction of 
the Ising variable s in H' allows us to perform Monte 
Carlo simulations. An (I + m)th order configuration 
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c = {s kl ,s k2 ,---,Sk n ; tk! , tk 2 j • ■ • i *fe„} is represented by 
the auxiliary spins s kl , Sfc 2 > ' " ' i s k n at the Keldysh times 
) tk 2 > " " j £fc n along the Keldysh contour, where the 
i(m) denotes the number of spins on the forward (back- 
ward) contour (see Fig. [1} and n = I + m. Its weight w c 
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FIG. 1: Illustration of the Keldysh contour for the CTQMC 
method. Arrows represent auxiliary Ising spins for a certain 
Monte Carlo configuration corresponding to the perturbation 
order I = 3 and m = 2 (n = 5). 



is then given as 
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where iV is an n x n matrix and its element consists of a 
2x2 matrix: 



U 
-(») 

with i, j = 1,2,- 
given by 



f («) - 



(11) 



= <Jfj exp (7s fci (T 2 ) , 
= &zGo(t ki7 t kj ), 
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where the times t' and t" correspond to the Keldysh times 
t' k and t k . The lesser and greater Green's functions are 
defined by a 2 x 2 matrix as 
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These Green's functions for the quantum dot system cou- 
pled to the normal and superconducting leads have been 
obtained by the standard technique , 11 ! 15 which are ex- 
plicitly represented in Appendix A. 

The sampling process must satisfy ergodicity and (as a 
sufficient condition) detailed balance. For ergodicity, it is 
enough to insert or remove the Ising variables with ran- 
dom orientations at random times to generate all possible 
configurations. To satisfy the detailed balance condition, 
we decompose the transition probability as 
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where pP ro P(p acc ) is the probability to propose (accept) 
the transition from the configuration i to the configura- 
tion j. Here, we consider the insertion and removal of the 
Ising spins as one step of the simulation process, which 
corresponds to a change of ±1 in the perturbation order. 
The probability of insertion/removal of an Ising spin is 
then given by 
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For this choice, the ratio of the acceptance probabilities 
becomes 
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where i) corresponding to a spin which is inserted 
on the forward (backward) contour. 

When the Metropolis algorithm is used to sample the 
configurations, we accept the transition from n to n ± 1 
with the probability 

i,q^jl. as) 

pP r °p(n± 1 ->• n)J 

In each Monte Carlo step, we can measure the Green 
function G(t,t'). By using Wick's theorem, the contri- 
bution of a certain configuration is given by 
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The expectation values of the particle number N(t) = 
J2a( n <?(t))> pairing correlations P(t) = (c^(i)c^(i)) , and 
double occupancy D(t) = (n^(t)n^(t)) at the quantum 
dot are calculated as 

N(t) = 2-i(G n (t,t))+i(G 22 (t,t)), (20) 

P(t) = i(G 12 (t,t)), (21) 

D(t) = l-N(t) + (detG(t,t)}. (22) 

We also measure the current from the quantum dot to 
the crth lead, which is given as 



I a = -21m Vkota (c kaa d a ) . 
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For convenient, we use the composite operator c aa = 
^2k VkairC kaa - and consider the following matrix as 
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The contribution for the matrix Aoa^t") of a certain 
configuration is given by 
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Then we obtain the measurement formula for the currents 
as 



I a (t) = -21m (A a (t,t)) n -(A a (t,t)) 
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This algorithm is essentially the same as the CTQMC 
method for the equilibrium stated and thereby it is 
straightforward to modify the codes to deal with the 
nonequilibrium system. Here, we comment on the dy- 
namical sign problem: the weight for a certain configu- 
ration is, in general, represented by the complex number 
[see cq. (fTU|) ]. As discussed in the previous works , 18 ' 19 
the dynamical sign problem becomes more serious in the 
simulations on longer contours and accurate measure- 
ments of physical quantities are restricted to a certain 
time t m ax- The introduction of the coupling strength Ts 
reduces the sign problem, which allows us to perform the 
simulations on longer contours. On the other hand, the 
bias voltage V little affects the dynamical sign problem. 
Therefore, performing Monte Carlo simulations with a 
fixed t max , we can equally treat the systems with differ- 
ent values of V, where the precision of the obtained re- 
sults little varies. This is contrast to the perturbative ap- 
proach, where more accurate results should be obtained 
in the vicinity of V = 0. 

In this study, we use the coupling constant of the nor- 
mal lead Tjv as the unit of energy, and fix the super- 
conducting gap as A/Tn — 0.5. In the following, we 
perform the CTQMC simulations to discuss nonequilib- 
rium behavior at zero temperature in the quantum dot 
coupled to the normal and superconducting leads. 



IV. RESULTS 

In this section, we discuss how the interaction quench 
affects the time evolution of the physical quantities. Fur- 
thermore, by extrapolating them in the t — > oo limit, we 



discuss steady-state properties of the quantum dot sys- 
tem. 

First, we consider the quantum dot system without 
the bias voltage. The time evolutions of the double oc- 
cupancy and pair correlation are shown in Figs. (a) 
and (c). In these figures, the quantities are shown on 
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FIG. 2: The double occupancy (upper panels) and pairing 
correlation (lower panels) in the system with T = 0, V = 0, 
A/rjv = 0.5 and td + U/2 — 0. The time evolutions are 
shown in (a) and (c), and the quantities at the steady state 
are shown in (b) and (d). 

the linear plot in the initial relaxation region (tT^ < 1) 
and on the logarithmic plot in the rest (£Tn > 1). When 
Ts = 0, the quantum dot is only coupled to the nor- 
mal lead and the system is reduced to the conventional 
Anderson impurity model, where our results reproduce 
the previous ones.— We find that the interaction quench 
decreases the double occupancy (tT^ < 1) and the sys- 
tem approaches the steady state (tV^ > 2). When the 
superconducting lead couples to the quantum dot, the 
double occupancy and pairing correlation for the initial 
state (t < 0) increase due to the proximity effect. As the 
interaction is turned on at t — 0, the double occupancy 
slightly decreases and the system quickly approaches the 
steady state, by comparing with the case Ts = 0, as seen 
in Fig. [2] (a). This implies that electron correlations 
become less important in the system. Although t max is 
finite, two quantities seem to converge around t = t max . 
This allows us to discuss the steady-state properties in 
the system. 

Regarding the physical quantities at t = t max as those 
for the steady state, we discuss the ground state prop- 
erties. The results are shown in Figs. [2] (b) and (d). 
When the quantum dot is only contacted to the nor- 



5 



mal lead (Ts = 0), the Coulomb interaction decreases 
the double occupancy and the Kondo-singlct dominant 
state is realized. As the coupling strength Ts increases, 
the double occupancy D and pair correlation P increase 
due to the proximity effect. In the large T$ region, the 
proximity-induced on-dot singlet-pairing dominant state 
is realized and electron correlations become irrelevant. In 
fact, the double occupancy approaches the value in the 
nonintcracting limit (D —> 1/4). Therefore, the crossover 
occurs between the Kondo-singlct and proximity-induced 
singlet-pairing dominant states.— 

When the bias voltage is applied, the current begins 
to flow between leads. In the steady state (t — > oo), 
the current / through the quantum dot is constant (I = 
In = —Is)- In the transient regime, the current from 
the dot to the crth lead is affected by the interaction U 
and hybridization T a , which results in different values at 
a certain time t. Therefore, we calculate both currents 
independently. The results for Ts/Tn = 1 are shown in 
Fig. [31 In the initial state (t < 0), the currents I a are 
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FIG. 3: The results are obtained by the CTQMC method 
in the quantum dot system with Fs/Tn = 1, V/Tn = 0.5, 
A/TV = 0.5 and T = when U/T N = 5.0 (circles) and 
U/Tn = 10.0 (triangles), (a) Solid (open) symbols rep- 
resent the time evolution of the current In (—Is)- (b) 
Solid (open) symbols represent the integrated total currents 

So, Jo 7a(i')<^' an d particle number N(t) at the quantum 
dot! 

given by the steady current through the noninteracting 
dot. The introduction of the interaction decreases the 
currents \In\ and \Ig\ differently, and oscillation behavior 
appears, as shown in Fig. [3] (a). Finally, Ijv ~ —Is and 



the system approaches the steady state. We note that the 
steady current can be extrapolated since the relaxation 
of oscillation behavior appears in the transient current. 

Here, we check the law of conservation of charge, which 
should be described as 



V I I a {t')dt' = N(t)-N(0). 
„ Jo 



(28) 



Fig. [3] (b) shows the integrated total currents from the 
quantum dot [left-hand side of eq. (|28p] and particle 
number at the quantum dot N(t). It is found that the dif- 
ference of two quantities is always constant, which means 
that the conservation law, cq. (|28|) , is satisfied within our 
numerical accuracy. This relation provides a good check 
for the numerical implementation. 
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FIG. 4: Current-voltage characteristics in the system with 
Ts/Tn = 1, A/Tjv = 0.5, and T = 0. Circles and triangles 
represent the CTQMC results for U/F N = 5.0 and 10.0 at 
t — tmax- Dashed lines represent the MPT results at the low 
temperature T = O.Oirjv.— 

Regarding the average of two currents at t = t max 
as the steady current [7 = (|Ijv| + |ls|)/2], we obtain 
the current-voltage characteristics in the systems with 
U/T N = 0.0,5.0, and 10.0, as shown in Fig. ffl As the 
error is defined as AI = \lN(tmax) — 1\, it is smaller than 
the size of symbols in the figure. When the bias voltage 
increases, the current monotonically increases, together 
with the kink structure around V = A. We also find 
that the increase of the Coulomb repulsion decreases the 
current. This implies that the interacting quantum dot 
can be regarded as the tunnel barrier on the interface. 
Although the maximum time in our simulations is lim- 
ited due to the dynamical sign problem (t max TN = 5.0 
for U/T N = 5.0 and i mox rjv = 2.0 for U/T N = 10.0), 
the CTQMC method reproduces reasonable results. In 
fact, in the weak coupling and small voltage region, our 
CTQMC data are in good agreement with the MPT re- 
sults at a very low temperature T/Tn = 0.01.— On the 
other hand, when V ~ A, the CTQMC data are away 
from the other. Since the precision of our data little 
depends on the bias voltage, this may suggest that the 
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MPT is not appropriate for quantitative discussions in 
this nonlinear response region V ~ A with U/Tn > 10. 

When the ratio Ts/rjv is away from unity, interesting 
behavior appears. The time evolutions of the currents for 
the systems with Ts/Fn = 0.2 and 5.0 are shown in Fig. 
In the Kondo-singlet dominant region (Fs/Fat = 0.2), 
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FIG. 5: (Color online) Time evolution of the currents for the 
systems with V/Y N = 0.5, A/Y N = 0.5, and T = 0. Open 
(closed) symbols represent the current from the superconduct- 
ing (normal) lead to the quantum dot when Ys/Yn = 0.2 (a) 
and Y S /Y N = 5.0 (b). 



tuates in the transient region iTjv < 1 and the system 
may not reach the steady state around t = t max . How- 
ever, its oscillation rapidly relaxes and the steady current 
is expected to be around two currents at t = t max . 

On the other hand, when the superconducting lead 
tightly couples to the quantum dot (Ts/Tn = 5.0), the 
currents are slightly changed by the interaction quench 
and the magnitudes of two currents approach each other 
around £Tjv ~ 0.5, in contrast to the Ts/Tn = 0.2 case. 
Nevertheless, we do not find the convergence in the cur- 
rent even around t = t max , as shown in Fig. [S] This may 
originate from the time evolution of the Kondo resonance 
induced by the interaction quench. It is characterized by 
an exponentially small energy scale and should be slow 
to be built up in the case Ts/^n > 1- In this case, 
the steady current, as which the transient currents at 
t = tmax are regarded in the paper, should be a lower 
bound of the correct value. Therefore, Monte Carlo sim- 
ulations on longer contours are necessary to obtain the 
steady current more precisely^ 

By performing similar calculations, we obtain the 
current-voltage characteristics of the systems with 
Ts/rAT = 0.2 and 5.0, as shown in Fig. [BJ In the nonin- 
teracting case (U = 0), the system should be regarded as 
the normal-superconducting junction with a simple tun- 
nel barrier. When Ts/^n is away from unity, the current 
rapidly increases around V ~ A. This behavior is well 
described by the conventional BTK theory^ The intro- 
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the interaction quench leads to a drastic change in the 
current In, in contrast to Is, as shown in Fig. [5] (a). 
When tTpf ~ 0.02, the current In changes its sign and a 
fairly large transient current flows against the bias volt- 
age around tTN ~ 0.3. This may be explained by the 
following. In the initial steady state, the particle number 
at the quantum dot is larger than the half filling. There- 
fore, the interaction quench tends to decrease the particle 
number at the quantum dot, which results in the decrease 
of both currents from the dots. In this case, there may be 
the relation AIs/AIn ~ r/v/rs around t ~ 0, which in- 
duces a large transient current between the quantum dot 
and the normal lead. As time progresses, the current In 
turns over again and the system approaches the steady 
state. In the case U/Tn = 10.0, the current largely fiuc- 



FIG. 6: Current-voltage characteristics in the system with 
rs/Tjv = 0.2 and 5.0 when A/Tjv = 0.5, e d + U/2 = and 
T = 0. Circles and triangles represent the CTQMC results 
for U/Ym = 5.0 and 10.0 at t — t max . Dashed lines represent 
the MPT results at the low temperature T = 0.0ir]v.— 

duction of the Coulomb interaction increases the current 
through the quantum dot tightly coupled to the super- 
conducting lead (Ts /Tn = 5.0). In the case, we could not 
obtain the steady current precisely, as discussed above. 
Nevertheless, our results are consistent with those ob- 
tained by the MPT,— which should be appropriate in 
this region. On the other hand, in the Kondo-singlet 
dominant region (Fs/Fat = 0.2), the Coulomb interac- 
tion decreases the current. In the MPT, the correlation 
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effects arc underestimated and the difference between two 
results becomes large in the strong coupling region. 

We finally discuss the effect of electron correlations in 
the system with a fixed voltage V/Ym = 0.5. The re- 
sults for U/Tn = 0.0,5.0, and 10.0 are shown in Fig. 
The increase in U monotonically decreases (increases) 




4 r s /r N 5 



FIG. 7: Currents functions of Ys/Yn in the system with 
V/Tjv = 0.5, A/Y N = 0.5, a d + U/2 = and T = 0. 

the current when r^/r^v < 1.5(r,g/rjv > 4). In the 
intermediate region (rs/IV ~ 2.5), nonmonotonic be- 
havior appears: the introduction of the interaction once 
increases the current in the proximity-induced on-dot 
pairing-dominant region. On the other hand, a further 
increase of the interaction decreases the current, where 
the Kondo singlet becomes dominant. The crossover be- 
havior in the nonlinear response region may be similar 
to that in the linear response region of the quantum dot 
system with A — > co^ where the zero bias conductance 
has a maximum around Ts ~ U/2. Therefore, we can 
say that the crossover between the Kondo-singlet and 
proximity-induced on-dot pairing dominant states affects 
the current- voltage characteristics under the high voltage 
beyond the superconducting gap V > A. 

In this paper, we have discussed the nonequilibrium 
phenomena of the quantum dot coupled to the normal 
and superconducting leads by means of the CTQMC 
method in the Keldysh formalism. It is an interesting 
problem how the local Coulomb interaction affects the 
Josephson current and multiple Andreev reflections in 
the quantum dot coupled to two superconducting leads, 
which is now under consideration. 



V. SUMMARY 

We have quantitatively studied the nonequilibrium 
phenomena through the quantum dot coupled to the 
normal and superconducting leads, extending the weak- 
coupling CTQMC method to treat the superconducting 
state in the Nambu formalism. We have confirmed that 
the obtained results are in good agreement with those 



obtained by the MPT in the weak coupling region and 
Ts/r^ > In the strong coupling region, we have 
clarified that the crossover between the Kondo singlet 
dominant and Cooper pairing singlet dominant regions 
affects the current- voltage characteristics. Transient cur- 
rents induced by the interaction quench have been dis- 
cussed. 
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VII. APPENDIX 

A. Calculation of noninteracting Green's functions 

Here, we explicitly show the impurity Green's func- 
tions in the noninteracting case U = 0. The Green's 
functions for the impurity (quantum dot) site are repre- 
sented in terms of the hybridization functions A (id) as, 



v</> 



(id ± irj)a - e d a z - A R ' A (oj) 



Gtf^M = G r (u)A</>(uj)G a (uj), 



where the hybridization functions are given as 

A R ' A (CU) : 



A«H 

Af(«) 



a<h 
a>h 



a=N,S 

~iY N a 
2iY N F{uj) 

2iY s Re \p(w)] f(u)M(u) 
-2iY N (a - F(u)] 



(29) 
(30) 

(31) 

(32) 
(33) 

(34) 

(35) 
(36) 
(37) 



A|(w) = -2ir s Re[ ( 8(id)][l-/(id)]M(o;), (38) 

where F(u) = diag [f(ui - fi N ), f(w + fi N )] , M(u>) = 
*o - ^a x and - 7 ^L Xf 6(\u\-A) - 
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\uj\). Using the Fourier transformations, we obtain the 
lesser and greater Green's functions eq. (|13l) . 

6</>(t) = / ^G< / >( W )e--. (39) 



where G^ d (i] = R,A,K) is the retarded, advanced, and 
Keldysh impurity Green's function for the Anderson 
model eq. flT]) and A% = A< + A> . Using the Fourier 
transformations as 



B. Calculation of Ao a (t,t') 

We obtain the expression of the matrix A 0a (t,t'), 
which is necessary to calculate the current from the quan- 
tum dot to the crth lead. We use the expression in pre- 
vious paper a 18 i 25 as 



2^ 



(41) 



aR aK 
M a 



^dd ^dd 

Gj d 



A R A K 



we obtain the matrix AQ a (t, t') in eq. ([24 



1 L. P. Kouwenhoven, D. G. Austing, and S. Tarucha, Rep. 
Prog. Phys. 64, 701 (2001). 

2 A. Eichler, M. Weiss, S. Oberholzer, C. Schonenberger, A. 
Levy Yeyati, J. C. Cuevas, and A. Martin-Rodero, Phys. 
Rev. Lett. 99, 126602 (2007). 

3 M. R. Graber, T. Nussbaumer, W. Belzig, and C. 
Schonenberger, Nanotechnology 15, S479 (2004). 

4 L. Hofstetter, S. Csonka, J. Nygard, and C. Schonenberger, 
Nature 461, 960 (2009). 

5 L. G. Herrmann, F. Portier, P. Roche, A. L. Yeyati, T. 
Kontos, and C. Strunk, Phys. Rev. Lett. 104, 026801 
(2010). 

6 R. S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, 
K. Shibata, K. Hirakawa, and S. Tarucha, Phys. Rev. Lett. 
104, 076805 (2010). 

7 R. S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, 
K. Shibata, K. Hirakawa, and S. Tarucha, Phys. Rev. B 81, 
121308 (2010). 

8 R. Fazio and R. Raimondi, Phys. Rev. Lett. 80, 2913 
(1998). 

9 P. Schwab and R. Raimondi, Phys. Rev. B 59, 1637 (1999). 
A. A. Clerk, V. Ambegaokar, and S. Hershfield, Phys. Rev. 
B 61, 3555 (2000). 

1 J. C. Cuevas, A. Levy Yeyati, and A. Martin-Rodero, Phys. 
Rev. B 63, 094515 (2001). 

2 T. Domahski, A. Donabidowicz, and K. I. Wysokihski, 
Phys. Rev. B 76, 104514 (2007). 



lx Y. Tanaka, N. Kawakami, and A. Oguri, J. Phys. Soc. Jpn. 

76, 074701 (2007). 
14 V. Koerting, B. M. Andersen, K. Flensberg, and J. Paaske, 

Phys. Rev. B 82, 245108 (2010). 

Y. Yamada, Y. Tanaka, and N. Kawakami, Phys. Rev. B 
84, 075484 (2011). 

16 J. Barahski and T. Domahski, Phys. Rev. B 84, 195424 
(2011). 

17 E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. 
Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011). 

18 P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B 79, 
035320 (2009). 

19 P. Werner, T. Oka, M. Eckstein, and A. J. Millis, Phys. 
Rev. B 81, 035108 (2010). 

20 E. Gull, P. Werner, O. Parcollet, and M. Troyer, Europhys. 
Lett. 82 57003 (2008). 

21 K. Satori, H. Shiba, O. Sakai, and Y. Shimizu, J. Phys. 
Soc. Jpn. 61, 3239 (1992). 

22 A. Koga and P. Werner, J. Phys. Soc. Jpn. 79, 064401 
(2010). 

23 G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Phys. 
Rev. B 25, 4515 (1982). 

24 A. F. Albuquerque et al., J. Magn. Magn. Mater. 310, 1187 
(2007). 

25 A.-P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 
50, 5528 (1994). 



